program stochastic_parameter_population

implicit none

interface
    subroutine random_normal(random_normal_1, random_normal_2)
    double precision :: random_normal_1, random_normal_2
    end subroutine
end interface

double precision :: theta, po, P, tw, ssf, tl, L, Sx, Sy, mean, sigma,random_normal_1, random_normal_2, meantw, sigmatw
double precision :: mean_ratio, sigma_ratio, mean_length, sigma_length, mean_prog, sigma_prog, mean_mig, sigma_mig
integer :: theta_int,P_int, i, tw_int, tl_int, L_int, bw_int, bl_int, Sx_int, Sy_int
character :: dummy1*1, filepath*1000
logical :: file_exists

call random_seed()

 write(*,*)"Enter the file path for the stochastic parameters output file"
 read(*,*) filepath 
  
inquire(file=trim(filepath)//"Stochastic parameters.txt", EXIST=file_exists)

if (file_exists)then;
open (unit=12,file=trim(filepath)//"Stochastic parameters.txt")
close(unit=12, STATUS='DELETE')
end if

theta=-1.0
do while (theta<0.0.or.theta>360.0)   
    call random_normal(random_normal_1, random_normal_2)
    
    write(*,*)"Enter mean and standard deviation for clinoform progradation angle"
        read(*,*) mean, sigma
    close(unit=1)
      
    if(random_normal_1<1)then;
        theta=mean+(sigma*random_normal_1)
    else if(random_normal_2<1)then;
        theta=mean+(sigma*random_normal_2)
    end if
    
    open(unit=2,file=trim(filepath)//"Stochastic clinoform progradation angle.txt")
    write(2,*)'!clinoform progradation angle'
    theta_int=theta
    write(2,*) theta_int
    close(unit=2)    
end do

po=-1.0
do while (po<0.0.or.po>1.0)     
    call random_normal(random_normal_1, random_normal_2)
    
    Write(*,*)"Enter mean and standard deviation for proximal boundary entry position"
        read(*,*) mean, sigma
    
    if(random_normal_1<1)then;
        po=mean+(sigma*random_normal_1)
    else if(random_normal_2<1)then;
        po=mean+(sigma*random_normal_2)
    end if

    open(unit=4,file=trim(filepath)//"Stochastic proximal boundary entry position.txt")
    write(4,*)'!Proximal boundary entry position'
    write(4,"(1F10.2)") po
    close(unit=4)   
end do
   
P=-1.0
do while (P<0.0)     
    call random_normal(random_normal_1, random_normal_2)
    
    write(*,*)"Enter mean and standard deviation of clinoform shape function exponent"
        read(*,*) mean, sigma
    
    if(random_normal_1<1)then;
        P=mean+(sigma*random_normal_1)
    else if(random_normal_2<1)then;
        P=mean+(sigma*random_normal_2)
    end if

    open(unit=6,file=trim(filepath)//"Stochastic clinoform shape function exponent.txt")
        write(6,*)'!clinoform shape function exponent'
            P_int=P
        write(6,*) P_int
    close(unit=6)    
end do

write(*,*)"Enter mean and standard deviation for the length of the top ellipse (tw)"
write(*,*)"in depostional strike"       
            read(*,*) meantw, sigmatw
            
write(*,*) "Enter mean and standard deviation for width/length ratio"
            read(*,*) mean_ratio, sigma_ratio
            
write(*,*)"Enter mean and standard deviation for clinoform length"
            read(*,*) mean_length, sigma_length
              
write(*,*)"Enter mean and standard deviation for clinoform progradation" 
            read(*,*) mean_prog, sigma_prog

write(*,*)"Enter mean and standard deviation for clinoform migration"
            read(*,*) mean_mig, sigma_mig
            
do i=1,2000
        tw=0
        ssf=-1
        L=0
        Sx=0
        Sy=0
    do while(tw<1)   
        call random_normal(random_normal_1, random_normal_2)
         
        if(random_normal_1<1)then;
            tw=meantw+(sigmatw*random_normal_1)
        else if(random_normal_2<1)then;
            tw=meantw+(sigmatw*random_normal_2)
        end if        
    end do
    
   do while(ssf<0)
        call random_normal(random_normal_1, random_normal_2)
    
        if(random_normal_1<1)then;
            ssf=mean_ratio+(sigma_ratio*random_normal_1)
        else if(random_normal_2<1)then;
            ssf=mean_ratio+(sigma_ratio*random_normal_2)
        end if       
    end do
    
    tl=tw/ssf
    
   do while(L<1)
        call random_normal(random_normal_1, random_normal_2)
    
        if(random_normal_1<1)then;
            L=mean_length+(sigma_length*random_normal_1)
        else if(random_normal_2<1)then;
            L=mean_length+(sigma_length*random_normal_2)
        end if       
    end do
    
      tw_int=tw
      tl_int=tl
      L_int=L
      bw_int=tw_int+L_int
      bl_int=tl_int+L_int
      
    do while(Sx<1)
        call random_normal(random_normal_1, random_normal_2)
   
        if(random_normal_1<1)then;
            Sx=mean_prog+(sigma_prog*random_normal_1)
        else if(random_normal_2<1)then;
            Sx=mean_prog+(sigma_prog*random_normal_2)
        end if       
    end do
    
    do while(Sy<1)
        call random_normal(random_normal_1, random_normal_2)
    
        if(random_normal_1<1)then;
            Sy=mean_mig+(sigma_mig*random_normal_1)
        else if(random_normal_2<1)then;
            Sy=mean_mig+(sigma_mig*random_normal_2)
        end if       
    end do
    
        Sx_int=Sx
        Sy_int=Sy
        
    open(unit=12,file=trim(filepath)//"Stochastic parameters.txt",ACCESS='APPEND')
    if (i==1)then;
        write(12,*)'!clinoform geometry and spacing parameters'
        write(12,*)'!bw bl  tw  tl  Sx  Sy'
    end if
        write(12,*) bw_int, bl_int, tw_int, tl_int, Sx_int, Sy_int
    close(unit=12)
end do
       
write(*,*) 'Press return to exit program'
read(*,*)

stop
end program 
!-----------------------------------------------------------------------------------------------------------------------------------------------
subroutine random_normal(random_normal_1, random_normal_2)

implicit none

double precision :: u, v, random_normal_1, random_normal_2, s
real :: pi
parameter (pi=3.14159265358979323846264338327950288419716939937510)

call random_seed()

s=2
random_normal_1=2
random_normal_2=2

do while(s>1) 
    call random_number(u)
    call random_number(v)
    s=(u**2)+(v**2)
end do

    random_normal_1=SQRT(-2*LOG(u))*(sin(2*pi*v))
    random_normal_2=SQRT(-2*LOG(u))*(cos(2*pi*v))   

return
end subroutine random_normal
!-----------------------------------------------------------------------------------------------------------------------------------------------
